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A hierarchical multiscale approach to model the magnetization dynamics of ferromagnetic ran¬ 
dom alloys is presented. First-principles calculations of the Heisenberg exchange integrals are linked 
to atomistic spin models based upon the stochastic Landau-Lifshitz-Gilbert (LLG) equation to 
calculate temperature-dependent parameters (e.g., effective exchange interactions, damping param¬ 
eters). These parameters are subsequently used in the Landau-Lifshitz-Bloch (LLB) model for 
multi-sublattice magnets to calculate numerically and analytically the ultrafast demagnetization 
times. The developed multiscale method is applied here to FeNi (permalloy) as well as to copper- 
doped FeNi alloys. We find that after an ultrafast heat pulse the Ni sublattice demagnetizes faster 
than the Fe sublattice for the here-studied FeNi-based alloys. 


I. INTRODUCTION 

Excitation of magnetic materials by powerful femtosec¬ 
ond laser pulses leads to magnetization dynamics on the 
timescale of exchange interactions. For elemental fer- 
romagnets the emerging dynamics can be p robed us¬ 
ing conventional magneto-optical methods 1 ®. For mag¬ 
nets composed of several distinct elements, such as fer- 
rimagnetic or ferromagnetic alloys, the individual spin 
dynamics of the different elements can be probed em¬ 
ploying ultrafast excitation in combination with the 
femtosecond-resolved x-ray magnetic circular dichroism 
(XMCD) techniqu^m. An astonishing example of such 
element-specific ultrafast magnetization dynamics was 
first measured on ferrimagnetic GdFeCo alloyJb There, 
it was observed that the rare-earth Gd sublattice demag¬ 
netizes in around 1.5 ps whereas the transition metal 
FeCo sublattice has a much shorter demagnetization time 
of 300 fs. Similar element-specific s pin dynamics was 
also observed in CoGd and CoTb alloyiP 7 T The element- 
selective technique allowed moreover to observe for the 
first time the element-specific dynamics of the so-called 
“all-optical switching” (AOSj® in GdFeCo alloys, find¬ 
ing that it unexpectedly proceeds through a transient- 
ferromagnetic-like state (TFLS) where the FeCo sublat¬ 
tice magnetization points in the same direc tion as that 
of the Gd sublattice before complete reversal 5 ®. Recent 
theoretical works supported the distinct demagnetiza¬ 
tion times observed experimentally®- 1 ® and their cru¬ 
cial role on the TFLS. AOS has been also demonstrated 
for other rare-earth transition-metal ferrimagnetic alloys 
as TbFtP, TbCcP, TbFeCcP, DyCcP, HoFeCcP, syn¬ 
thetic ferrimagnet!®3®] an d very recently in the hard- 
magnetic ferromagnet FePt 19 . 

Although the full theoretical explanation of the 
therm ally driven AOS process is still a topic of 
debateM3ME3 the distinct demagnetization rates of 


each of the constituting elements has been suggested as 
the main driving mechanism for the A OS observed on 
antiferromagnetically coupled alloy * -* — These findings 
have highlighted the question how ultrafast demagneti¬ 
zation would proceed in ferromagnetically coupled two- 
sublattice materials such as permalloy (Py). Unlike rare- 
earth transition-metal alloys which consists of two intrin¬ 
sically different metals, Py is composed of Fe (20 %) and 
Ni (80 %) which have a rather similar magnetic nature, 
due to a partially filled 3d shell. Thus, it is a priori not 
clear if their spin dynamics should be the same or differ¬ 
ent. 

Recent measurements have addressed this question. 
Using extreme ultraviolet pulses from high-harmonic 
generation sources Mathias et al ^ 1 probed element- 
specifically the ultrafast demagnetization in Py and ob¬ 
tained the same demagnetization rates for each element, 
Fe and Ni, but with a 10 to 70 fs delay between them. 

From a theoretical viewpoint an important question 
is which materials parameter are defining for the ultra¬ 
fast demagneti zation . Thus far, different criteria have 
been suggeste d 25 * 26 ! For single-element ferromagnets, 
Kazantseva et estimated, based on phenomenolog¬ 
ical arguments, that the timescale for the demagnetiza¬ 
tion processes is limited by ss p/(2A^k-QT pn \ se ). 

Here, 7 d ema g depends not only on the elemental atomic 
magnetic moment, /i, but also on the electron tempera¬ 
ture, Tp U i se , and on the damping constant A. Assuming 
that the damping constants A and gyromagnetic ratios 7 
are equal for Fe and Ni the demagnetization time would 
therefore only vary due to the different magnetic mo¬ 
ments of the constituting elements. In that case, the 
demagnetization time of Fe is larger than the one for Ni 
(since p Fe > p Nl , see Table [i] below). 

A similar criterion (as in Ref. [23 for single-element fer¬ 
romagnets) has been suggested by Koopmans et alW^ on 
the basis of the ratio between the magnetic moment and 
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the Curie temperature, h/Tq. Since for ferromagnetic al¬ 
loys each element has the same Curie temperature, this 
criterion would lead to the same conclusions as Kazant¬ 
seva et al.', the different atomic magnetic moments of Fe 
and Ni are responsible for the different demagnetization 
times. Furthermore, Atxitia et al. 10 have theoretically 
estimated the demagnetization times in GdFeCo alloys 
proposing that the demagnetization times scale with the 
ratio of the magnetic moment to the exchange energy of 
each element and a similar relation is expected for ferro¬ 
magnetic alloys. The demagnetization times of Fe and Ni 
in Py were also theoretically investigated by Schellekens 
and Koopmans in Ref. QT] where a modified microscopic 
three temperature model (M3TM)®was used. Thereby, 
they obtained a perfect agreement with experimental re¬ 
sults of Mathias et «/.j® but only assuming an at least 
4 times larger damping constant for Fe. However, this 
work does not provide a simple general criterion, valid 
for other ferromagnetic alloys. 

We have developed a hierarchical multiscale approach 
(cf. Ref. E7j) to investigate the element-specific spin dy¬ 
namics of ferromagnetic alloys and to obtain a deeper 
insight into the underlying mechanisms. First, we con¬ 
struct and parametrize a model spin Hamiltonian for 
FeNi alloys on the basis of first-principles calculations 
[Sec. IIA . This model spin Hamiltonian in combina¬ 
tion with extensive numerical atomistic spin dynamics 
simulations based on the stochastic LLG equation are 


used to calculate the equilibrium properties [Sec. IIB 


as well as the demagnetization process after the appli¬ 
cation of a step heat pulse. The second step of the pre¬ 
sented multiscale model links the atomistic spin model to 
the macroscopic two-sublattices Landau-Lifshitz-Bloch 
(LLB) equati on of motion recently derived by Atxitia 
et alW - 1 [Sec. Ill. The analytical LLB approach allows 


for efficient simulations, and most importantly, provides 
insight in the element-specific demagnetization times of 
FeNi alloys. 


II. FROM FIRST PRINCIPLES TO ATOMISTIC 
SPIN MODEL 

A. Building the spin Hamiltonian 

To start with, we construct an atomistic, classical spin 
Hamiltonian T~L on the basis of first-principles calcula¬ 
tions. In particular, we consider three relevant alloys: 
Fe 5 oNi 5 o, Fe 20 Ni 8 o (Py) and Py 60 Cu 4 o. The first two al¬ 
loys will allow us to assess the influence of the Fe and Ni 
composition, while the last two alloys will permit us to 
study the effect of the inclusion of non-magnetic impuri¬ 
ties on the demagnetization times. This was motivated 
by the work of Mathias et aZ.® who studied the influence 
of Cu doping on the Fe and Ni demagnetization times in 
an Py 60 Cu 40 alloy. 

To obtain the spin Hamiltonian we have employed spin- 
density functional theory calculations to map the behav¬ 


ior of the magnetic material onto an effective Heisenberg 
Hamiltonian, which can be achieved in various ways 29 30 . 
Here we use the two-step approach suggested by Licht¬ 
enstein et The first step represents the calculation 

of the self-consistent electronic structure for a collinear 
spin structure at zero temperature. In the second step, 
exchange parameters of an effective classical Heisenberg 
Hamiltonian are determined using the one-electron Green 
functions. This method has been rather successful in ex¬ 
plaining magnetic thermodynamic properties of a broad 
class of magnetic materials®® 

The self-consistent electronic structure was calcu¬ 
lated using the tight-binding linear muffin-tin orbital 
(TB-LMTO) approach 32 within the local spin-density 
approximation 3 ^ to the density functional theory. 

Importantly, the materials we investigate here are al¬ 
loys. Hence, it is assumed that atoms are distributed 
randomly on the host fee lattice. The effect of disor¬ 
der was described by the coherent-potential approxima¬ 
tion (CPA)®. The same radii for constituent atoms were 
used in the TB-LMTO-CPA calculations. We have used 
around a million ^-points in the full Brillouin zone to 
resolve accurately energy dispersions close to the Fermi 
level. 

The calculations of the Heisenberg exchange constants 
Jij in ferromagnets can be performed with a reason¬ 
able nu meric al effort by employing the magnetic force 
theorerrPsI 3 ^ It allows to express the infinitesimal 
changes of the total energy using changes in one-particle 
eigenvalues due to non-self-consistent changes of the ef¬ 
fective one-electron potential accompanying the infinites¬ 
imal rotations of spin quantization axes, i.e., without any 
additional self-consistent calculations besides that for the 
collinear ground state. The resulting pair exchange in¬ 
teractions are given by 

Jij = -Im [ d IE /dr [ dr'B ex (r) GlB e ^ (r') G l _, (1) 
^ J —oo J Qi J j 

with G\ = G^(r,r',E+) and G l _ = G^ (r',r,E~). E F 
denotes the Fermi level and fb the i-th atomic cell, 
a =t,i is the spin index, E + = limc^o-E + ia, G a 
are spin-dependent one-electron retarded Green func¬ 
tions, and Hex is the magnetic field from the exchange- 
correlation potential. The validity of this approxima¬ 
tion has been examined more quantitatively in several 
studies. 3 ^® The ab initio calculated distance-dependent 
exchange constants for the Fe 2 oNi§o alloy, i.e., the ex¬ 
change within the Fe sublattice (Fe-Fe), the Ni sublat¬ 
tice (Ni-Ni) as well as between the Fe and Ni sublattices 
(Fe-Ni), are shown in Fig. [lj The calculated magnetic 
moments for all three alloys considered here are given in 
Table H 

In our hierarchical multiscale approach, these com¬ 
puted material parameters (the exchange constant ma¬ 
trix as well as the magnetic moments) are now used as 
material parameters for our numerical simulations based 
on an atomistic Heisenberg spin Hamiltonian. We con¬ 
sider thereto classical spins Sf = A 4 ! /with e randomly 
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TABLE I. Ab initio calculated magnetic moments and experimental lattice constants A used in the atomistic Langevin spin 
dynamics simulations. Effective exchange parameters calculated from ab initio calculations, Jg S = JA Jqj, where the sum is 

here over all neighbors j. Curie temperatures as calculated from the atomistic simulations, Tq LG , and the experimental value, 

rj- iexp 


alloy 

M Fe ia m A 

[Mb] [mb] [nm] 

rNi-Ni 

J 0 

[J xHT 21 ] 

TFe —Fe 
J 0 

[Jxl0~ 21 ] 

rFe-Ni 

J 0 

[J x 10~ 21 ] 

rpLLG 

1 C 

[K] 

rj-iex. p 

[K] 

Py 

2.637 0.628 0.355(P 

6.2419 

32.3162 

26.3654 

650 

85(P 

NisoFeso 

2.470 0.730 0.358^ 

6.6265 

25.3789 

25.0656 

850 


PyeoCruo 2.645 0.429 0.3550 

2.6623 

56.2789 

22.6442 

340 

40fP 



FIG. 1. (Color online) Ab initio calculated exchange constants 
Jij for the Fe 2 oNi 8 o alloy of the distance i'ij between atoms 
i and j. Results are given for the three different possible 
sublattice interactions (JFe-Fe, Jm—au and JFe-Ni). Note 
our hyperbolic scaling. In our atomistic spin simulations the 
exchange constants are taken into account up to a distance rij 
(cutoff) where they are finally small enough to be neglected. 


representing iron (/if = /tf e ) or nickel magnetic moments 
(/if = /if 1 ) on the fee sublattice. For the Cu-doped 
Py 6 oCu 4 o alloy the calculated magnetic moments on Cu 
vanish, i.e. /xf u = 0. 

The spin Hamiltonian for unit vectors, Sf, representing 
the normalized magnetic moments of the *-th atom on 
either the Fe or Ni sublattice reads 

w = -E(T s f' s ' (2) 

ij 

Mo/jf/4 3(Sf ■ ejj)(ejj ■ S'-) — Sf ■ S^ \ 

Stt r?- J’ 

The first sum represents the exchange energy of mag¬ 
netic moments, either on Ni or on Fe sites, distributed 
randomly with the required concentrations. The ex¬ 
change interaction matrices (corresponding to Jm-Nn 
t/pe—Ni, or JNi-Ni) are those from the ab initio calcula¬ 


tions (as shown for Py in Fig.[lJ. These have been taken 
into account up to a distance of six unit cells (cutoff also 
shown in Fig. [l]) until they are finally small enough to 
be neglected. The second sum describes the magnetic 
dipole-dipole coupling. 

Note, that the exchange interaction given by the ma¬ 
trices Jij is incorporated in our atomistic spin dynamics 
simulations via the Fast Fourier transformation method 
(see Ref.SZIfor more details). As a side effect, we are able 
to calculate the dipolar interaction without any addi¬ 
tional computational effort so that we take them into ac¬ 
count although they will not influence our results much. 

Since we are interested in thermal properties we 
use Langevin dynamics, i.e. numerical solutions of the 
stochastic LLG equation of motion 

C 1 + ( y K Sf = -Sf x [Hi + Af (Sf x HO], ( 3 ) 

li 

with the gyromagnetic ratio 7 ?, and a dimensionless 
Gilbert damping constant Af that describes the coupling 
to the heat-bath and corresponding either to Fe or to Ni. 
Thermal fluctuations are included as an additional noise 
term (i in the internal fields H; = — + C,i{t) with 

(Ci(t)) = 0, (C»T,(0)OeW) = B / j ^M(«). (4) 

(i 

where i,j denotes lattice sites occupied either by Fe or 
Ni and 77, 9 are Cartesian components. All algorithms we 
use are described in detail in Ref. FCT1 


B. Equilibrium properties: element-specific 
magnetization 

First, we investigate the element-specific zero-field 
equilibrium magnetizations for Fe and Ni sublattices. 
Those magnetizations are calculated as the spatial and 
time average of the sum of local magnetic moments, 
m e = (S e ) with e representing either Fe or Ni. For 
our numerical studies, we assume identical damping con¬ 
stants (A = Af) as well as gyromagnetic ratios (7 = 7f 
= 1.76 • 10 11 (Ts) -1 ) for both, Fe or Ni. We perform 
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FIG. 2. (Color online) Element-specific zero-field equilibrium 
magnetizations m e of either Fe or Ni as a function of tem¬ 
perature calculated by a rescaled mean-field approximation 
(MFA) (lines) and by the atomistic spin dynamics simulation 
(open symbols). In the MFA the exchange parameters are 
renormalized by equalizing the Curie temperatures Tc com¬ 
puted with atomistic simulations with those obtained from 
the rescaled MFA. System size 128 x 128 x 128, damping 
parameter A = 1.0. 


our Langevin spin dynamics simulations for two differ¬ 
ent FeNi alloys, namely FesoNiso and Py, as well as for 
permalloy diluted with copper, Py 6 oCu 4 o- All material 
parameters used in our simulations are given in Table [I] 


The temperature dependence of the normalized 
element-specific magnetizations m e are shown in Fig. 
[2] The calculated values of the Curie temperatures are 
given in Table |T] together with known experimental val¬ 
ues. Both, the numerical and experimental values, are 
in good agreement. The element-specific magnetizations 
as well as the total magnetization (not shown in Fig. [2]) 
of the alloys share the same Curie temperature while in 
the temperature range below the Curie temperature their 
temperature dependence is different for the two sublat¬ 
tices; the normalized magnetization of Ni is lower than 
that of Fe. 


The element-specific magnetizations calculated within 
the framework of a rescaled mean-field approximation 
(MFA) are shown as well. This approach will be dis¬ 
cussed in detail in Sec. IIIII below where these curves serve 
as material parameters for the simulations based on the 
LLB equation of motion also introduced in the next sec¬ 
tion. 


III. FROM ATOMISTIC SPIN MODEL TO 
MACROSCOPIC MODEL 


A. Two-sublattices Landau-Lifshitz-Bloch equation 


Within the hierarchical multiscale approach, the 
macroscopic (micromagnetic) equation of motion valid at 
elevated temperatures is the LLB equation 27 ! Initially, 
the macroscopic LLB equation of motion was derived by 
Garanin for single-species ferromagnets only. Garanin 
first calculated the Fokker-Planck equation for a single 
spin coupled to a heat-bath, thereafter a non-equilibrium 
distribution function for the thermal averaged spin polar¬ 
ization was assumed to drive the non-equilibrium dynam¬ 
ics. Second, the exchange interactions between atomic 
spins were introduced using the mean field approximation 
(MFA) with respect to the spin-spin interactions. This 
last step reduces to the replacement of the ferromagnetic 
spin Hamiltonian % with the MFA Hamiltonian "H MFA in 
the single (macro)spin solution. 

The LLB formalism was recently broadened to de¬ 
scribe the distinct dynamics of two-sublattices mag¬ 
nets, both antiferromagnetically or ferromagnetically 
coupled®!. The derivation of such equations follows sim¬ 
ilar steps as for the ferromagnetic LLB version but con¬ 
sidering sublattice specific spin-spin exchange interac¬ 
tions and MFA exchange fields, (H^ PA ) conf . For the ex¬ 
change field the random lattice model is used by gener¬ 
ating the random average with respect to disorder con¬ 
figurations (... )° ont . The corresponding set of coupled 
LLB equations for each sublattice reduced magnetiza¬ 
tion m e = (S e ) = IVL/ilLf, where M s e is the saturation 
magnetization at 0 K, has the form 


riff = 7 e K x (H^ fa 

— Ff, | 1 — m mr 


- n 


m x nr x m , 


(m e ) 2 


(m e ) 2 


m . 


(5) 


Here, mg = £(£o)§^ is the transient (dynamical) magne¬ 
tization to which the non-equilibrium magnetization m e 
tends to relax, and where £g = Ur (Hmfa)™' is the ther¬ 
mal reduced field, £g = |£g|, and C (£) = coth (£) — l/£ is 
the Langevin function and £'(£) = d£(£)/d£;. The par¬ 
allel (rjj) and perpendicular (r^) relaxation rates in Eq. 
© are given by 


FS=A( 


1M 

a £'(£§} 


and 




( 6 ) 


A^j = 2ksT^ \ e /is the characteristic diffusion relax¬ 
ation rate. The damping parameters A c have the same 
origin as those used in the atomistic simulations. 

The first and the second terms on the right-hand side 
of Eq. ([5]) describe the transverse motion of the mag¬ 
netization. These dynamics are much slower than the 
longitudinal magnetization dynamics given by the third 
term in this equation. Therefore, in the following we will 
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FIG. 3. (Color online) Schematics of the magnetic unit cell 
used in the mean-field approximation for the FeNi alloys. The 
unit cell shown by the box contains two spins, one Fe and one 
Ni. The only interaction among spins located at the same unit 
cell r is defined by J^~ Fe . The self-interactions are neglected, 
jNi-Ni( r , _ jFc- c (r,r) = 0. The rest of the interactions 
are among spins located in neighboring unit cells r and r' . 


neglect the transverse components (in Eq. ([5])) and keep 
only the longitudinal one, 

m e = —rf| (to 6 - TOq) . (7) 

In spite of the fact that the form of Eq. (J7j) is similar 
to the well known Bloch equation, the quantity too = 
?7 io (to 6 , to^) (with S the 2-nd type of element) is not 
the equilibrium magnetization but changes dynamically 
through the dependence of the effective field (H^ FA ) oonf 
on both sublattice magnetizations. Moreover, the rate 
parameter Tf, = rjj (toq,TOq) contains highly non-linear 

terms in toJ) and TOq. 

Therefore, the analytical solution of Eq. ([7| and thus 
a deeper physical interpretation of the relaxation rates is 
difficult without any further approximations. However, 
Eq. (J7| can be easily solved numerically with the aim to 
directly compare the solutions to those of the atomistic 
spin simulations. This is discussed in more detail in the 
next subsections. 


B. FFom atomistic spin model to 
Landau-Lifshitz-Bloch equation 

Next, to solve Eq. ([5]) or Eq. ([7]), one needs to calculate 
(H^ FA ) conf for the here-considered FeNi alloys. An ade¬ 
quate definition of such a field will allow us to directly 
compare the magnetization dynamics from our atomistic 
spin simulation with the LLB macroscopic approach. 

However, a quantitative comparison between both a 
standard MFA and atomistic spin model calculations of 
the equilibrium properties is usually not possible. This is 
due to the fact that the Curie temperature gained with 
the MFA approach is overestimated due to the inher¬ 
ent poor approximation of the spin-spin correlations. Al¬ 
though, rescaling the exchange parameters conveniently 
in such a way that the Curie temperature calculated 
with the MFA approach agrees with atomistic simula¬ 
tions leads to a good agreement of both methods. Hence, 


we first present the standard MFA for disordered two- 
sublattices magnets, thereafter, we will deal with the 
rescaling of the exchange parameters. 

The MFA Hamiltonian of the full spin Hamiltonian for 
FeNi alloys (see Eq. ([2]) introduced in Sec. 0 can be 
written as 

Hmfa = H 00 -// e E H- FA .Sr-M Ni E H MPA-Sr, (8) 

i i 

where the dipolar interaction is neglected. The mean 
field acting on each site i can be separated in two contri¬ 
butions; a) the contribution from neighbors of the same 
type j e an d b) those of the other type j s , 

// (HL FA r nf = e 4< s 4 + E 4< s 4> (9) 

ej‘ ej 6 

where sums run over the nearest neighbours. When the 
homogenous magnetization approximation is applied (i.e. 
(SjPe) = m Fe and (SjNi) = m Ni for all sites) one can de¬ 
fine ,Jq 4 anc ^ 4 = X4 4- A sketch of the 

exchange interaction within the present MFA model is 
presented in Fig. [3] The impurity model is mapped to a 
regular spin lattice where the unit cell (orange box) con¬ 
tains the two spin species, Fe and Ni, and the exchange 
interactions among them are weighted in terms of the 
concentration of each species. 

The equilibrium magnetization of each sublattice to 6 
can be obtained via the self-consistent solution of the 
Curie-Weiss equations = £(4r (Hi/fa) 00 ”*)- 

Fig. [2] shows good agreement of the calculated m%{T) 
using the MFA and the atomistic spin model for the 
three system studied in the present work. The exchange 
interactions are rescaled as J(f Mfa ~ (1.65/2) Jq 5 , for 
FesoNiso and Py. For Py 6 oCu 4 o it is in agreement with 
4 mfa = (1-78/2) Jg* 5 . Here, the atomistic calculations is 
not as accurate for intermediate temperatures as for the 
other two alloys. This could be because of the increased 
complexity introduced by the inclusion of Cu impurities 
which cannot be fully described by the MFA. 

C. De- and remagnetization due to a heat pulse 

In the following, we study the reaction of the element- 
specific magnetization to a temperature step in Py as well 
as in Py diluted with Cu. In the first part of the temper¬ 
ature step the system is heated up to T = 0.8 Tq and in 
the second part it is cooled down to T pu i se = 0.5 Tq. The 
heat pulse roughly mimics the effect of heating due to a 
short laser pulse. The first part of the temperature step 
triggers the demagnetization while the second one trig¬ 
gers the remagnetization process. We perform atomistic 
as well as LLB simulation of the de- and remagnetization 
of the two sublattices after the application of a step heat 
pulse of 500-fs duration. 

The reaction of the Fe and Ni sublattice magnetiza¬ 
tions is shown in Fig. [4] While the temperature step 
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FIG. 4. (Color online) Calculated z-component of the nor¬ 
malized element-specific magnetization mj vs. time for Py 
(top panel) and Py 6 oCu 4 o (bottom panel). In both cases the 
quenching of the element-specific magnetizations for Fe and 
Ni due to a temperature step of T pu i se = 0.8 Tc are shown, 
computed with atomistic Langevin spin dynamics (open sym¬ 
bols) as well as LLB simulations (lines). System size 64 x 64 
x 64, damping parameter A = 0.02. 


is switched on, the two sublattices relax to the corre¬ 
sponding equilibrium value of the sublattice magnetiza¬ 
tions m e (Tp U i se ). Note, that these equilibrium values 
are different for the two sublattices in agreement with 
the temperature-dependent equilibrium element-specific 
magnetizations shown in Fig. [2j 

Because of that, the different demagnetization time 
scales are not well distinguishable in Fig. [4j Thus, 
we use the normalized magnetization, m porm = (m e — 
m min)/( m (t=o) — m min) °f the sublattices, rather than 
m e to directly compare the demagnetization times. The 
demagnetization time after excitation with a tempera¬ 
ture pulse Master for Ni than for Fe (Fig. [5] (top panel)) 
for the first 200 fs, while one can see that for times larger 
than 200 fs both elements demagnetize at the same rate 
(Fig. [5] (bottom panel)). Experiments on Py suggest that 
the time shift between distinct and similar demagnetiza¬ 
tion rates in Py is of around 10-70 femtoseconds®! 


D. Understanding relaxation times within the 
Landau-Lifshitz-Bloch formalism 

The relaxation rates of the Fe and Ni sublattices 
can be understood by discussing the linearized form of 
Eq. 0. Here, the expansion of Tjj and m g around 
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FIG. 5. (Color online) Top panel: Normalized magnetiza¬ 
tion dynamics of Fe and Ni sublattices after the application 
of a heat pulse T = 0.8 Tc as computed with the atomistic 
spin model. The ratio between the Fe and Ni demagneti¬ 
zation times is 1.8. The intersection of the linear fit to the 
abscissa gives the relaxation time for each sublattice. Bot¬ 
tom panel: plot of the unnormalized magnetization dynamics 
which shows that after the first 0.2 picosecond the element- 
specific demagnetization proceeds at the same rate. 



FIG. 6. (Color online) Relaxation times of the dynamical 
system obtained by the LLB equation as a function of tem¬ 
perature. Inset: The ratio between the relaxation times. 


their equilibrium values is considered®! and leads 
to <9(Am )/dt = ^4||Am with Am = (Am 6 , Am s ) and 

m e(<5) _ jTjd 5 ) _|_ A m £ ( s \ Furthermore, the characteristic 
matrix A\\ drives the dynamics of this linearized equation 
and has the form 


A\\ = 


- 7 e a|j/A“ \ 

J 0 5 7/ -7^f/A 55 ) ’ 


with 


( 10 ) 


; Jp m e e 
jA m s e 


and 


A“ = 


Xy 


1 + tM 


( 11 ) 
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where \\\ are the longitudinal susceptibilities which can 
be evaluated in the MFA approximation as 

_ Jt*v, s C s C* + ^C<(k B T - J*£ s ) 

X|1 {k B T- J*C s )(k B T- J 0 e £ £ )- 1 1 


with C e = £'(£!) and C 5 = d (£f). We note that the 
longitudinal susceptibility in Eq. (121 depends on the ex¬ 
change parameter (Curie temperature) and the atomic 
magnetic moments of both sublattices. 

Next, the longitudinal damping parameter in Eq. (10) 
is defined as a e = {2k B TX L m e e ) /ii e H^ ex , where is 

the average exchange field for the sublattice e at equi¬ 
librium, defined by the MFA expression The longi¬ 
tudinal fluctuations are defined by the exchange energy, 
according to the expression above. However, the longitu¬ 
dinal relaxation time is not simply inversely proportional 
to the damping parameter. Instead the relaxation pa¬ 
rameters in Eq. (10) do also depend on the longitudinal 
susceptibilities which give the main contribution to their 
temperature dependence. 

It is important to note that the matrix elements in 
Eq. (101 are temperature as well as (sublattice) material 
parameter dependent. The general solution of the char¬ 
acteristic equation, |A|| — T ± I\ = 0, gives two different 


eigenvalues, = 1 /t±, corresponding to the eigenvec¬ 
tors v±. Here, I is the unit matrix. The computed tem¬ 
perature dependence of the relaxation times t± is pre¬ 
sented in Fig. [g] More interestingly, we observe that the 
ratio between relaxation times r + /r_ [inset Fig. [fi] is al¬ 
most constant for temperature below 0.5 Tq and it has a 
value of 1.8 which compares well with atomistic simula¬ 
tions [Fig. [5]. At elevated temperatures, one relaxation 
time r + will dominate the magnetization dynamics of 
both sublattices. 

In Fig. [7](a) we present the temperature dependence of 
the longitudinal damping parameters and in Fig. mb) 
the temperature dependence of the parameters r e ° = 
ajj/A e<5 . These parameters define the element-specific 
longitudinal dynamics. In Figs. [7] (c) and (d) the temper¬ 
ature dependent a[j/ay and A ee /A ss are shown. It can 
be seen that at least in the range of low temperatures the 
magnetization dynamics is mainly defined by r ee T 65 . 

The general solution of the linearized LLB system for 
the two sublattices can be written as 


Am Fo (t) = A Fe exp (— t/r + ) + B Fc exp (—f/r_) 
Am Nl (t) = A Nl exp (—t/r+) + B Nl exp (—i/r_), (13) 

where the coefficients A Fe ( N d and H Fe ( N d will depend 
of the eigenvectors v± and the initial magnetic state, 
Am Fe (0) and Af7i Nl (0). For instance 


A Fe 


Am Fe (0) 


Am Ni (0) 
Am Fe (0) X + 


— X+ 



(14) 


where x + = u Fe /u^' and x_ = u Fe /u^‘, is the ratio be¬ 
tween he eigenvector components. The other coefficients 
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FIG. 7. (Color online) (a) Temperature dependence of the in¬ 
dividual longitudinal damping parameters for Fe and Ni. (b) 
Matrix elements of the dynamical system defining the magne¬ 
tization dynamics, (c) Ratio between the individual damping 
parameters, (d) Ratio between the “effective” susceptibilities 
A“ and the actual susceptibilities \'[| • 


are calculated similarly. This complexity prohibits a gen¬ 
eral analysis of the results. Thus, although the general 
solution is clearly a bi-exponential decay, one can wonder 
when the one exponential decay approximation will give 
a good estimate for the individual relaxation dynamics. 

Two interesting scenarios exist: First, the relaxation 
times t + and r_ could have very different time scales 
and thus one can separate the solution on short and long 
time scales, defined by r_ and r+, respectively. This 
is an interesting scenario for ultrafast magnetization dy¬ 
namics where only the fast time scale will be relevant. 
Fig. § shows the ratio t + /t_ and we can observe that 
the scenario r+/r_ 1 only happens for temperatures 

approaching Tq . As we have seen in the atomistic simula¬ 
tions, after an initial distinct quenching of each sublattice 
magnetization, both sublattice demagnetize at the same 
rate but slower than the initial rates (see Fig. [ 5 ]). 

The second scenario occurs when A Fe « Am Fc (0) and 
B m « Am Ni (0), even if r + and r_ are of the same or¬ 
der. This happens, for example, either when the coupling 
between sublattices is very weak, or at relatively low tem¬ 
peratures, see Fig. [G] In this case the system can be con¬ 
sidered as two uncoupled ferromagnets (although with 
renormalized parameters), meaning that the matrix in 
Eq. (10) defining the dynamics is almost diagonal. Thus, 
we can approximately associate each eigenvalue of Eq. 
(101 to each sublattice, r_ = r Nl and r + = r Fe . The in¬ 
set in Fig. [6] shows the ratio r+/r_ for the whole range of 
temperatures. At low-to-intermediate temperatures we 
find that r + /r_ « 1.8. This is in good agreement with 

























TABLE II. Theoretical results: ab initio calculated ratio be¬ 
tween the mean exchange interaction at T = 0 K, the ratio 
between atomic magnetic moments and the quotient of these 
ratios. Results of simulations: atomistic spin model calcu¬ 
lated ratio between k exponents and relaxation times. The 
ratio between the magnetic atomic moments and the expo¬ 
nents k is predicted in the main text to give the ratio between 
relaxation times. 


alloy 

theoretical 

simulations 

?Fe 

J 0 

jNi 

J 0 

e Pe 

M Ni 

Fe /Ni 

M J 0 

« p<> 

r Fe 

7NT 

M Fe K Pe 
,*Ni re Ni 

Fe 5 oNi 5 o 

1.592 

3.38 

2.12 

1.492 

2.10 

2.25 

Py 

2.685 

4.198 

1.563 

2.3 

1.8 

1.8 

Py6oCu4o 

4.412 

6.17 

1.398 

2.95 

2.1 

2.05 


atomistic simulations, see Fig. |5ja), and it clearly shows 
that the relaxation times ratio is not defined by the ratio 
between atomic magnetic moments, /r Fo //i Nl ss 4. 

In the case that the longitudinal relaxation rates are 
defined by the diagonal elements of the matrix (10 ) and 
T is not close to Tq the longitudinal relaxation time can 
be estimated as 


1 


2y e A e m%Hl 


(15) 


Thus the ratio between the relaxation rates of Ni and Fe 
(for the same gyromagnetic ratio value, the same cou¬ 
pling parameter and not too close to Tq) is defined by 


r Ni _ /A Fe /z Ni \ J Fe m Fe 


(16) 


We recall that J^ml = Jqto| + Jo S mi is the average ex¬ 
change energy for the sublattice e at equilibrium. Thus, 
the interpretation of the ratio of the relaxation times is 
straightforward. The low temperature value of the ra¬ 
tio Jq°/ Jq 1 is presented in Table [n] for the three alloys 
studied here. The second column presents the ratio be¬ 
tween atomic magnetic moments, and the third column 
the estimated ratio between relaxation times under the 
assumption of equal damping parameter at each sublat¬ 
tice. 

The estimated ratios for relaxation times are in rather 
good agreement with the atomistic simulations (fifth col¬ 
umn) for Fe 5 oNi 50 and Py, however for Py 60 Cu 4 g the es¬ 
timation is not that good. We have to remember that the 
MFA re-scaling of the exchange parameters did not give 
a completely satisfactory result for the shape of m(T) in 
this alloy (see Fig. 2ja)). Thus, since the re-scaled ex¬ 
change parameter does not work completely well at the 
low-to-intermediate temperature interval, we further in¬ 
vestigate this case (Py 6 oCu 4 o) by relating the obtained 
relation in Eq. ( fTtij ) for the ratio T Nl /r Fe to the slopes of 
the curves m{T). 


This can be easily done by using the linear decrease of 
magnetization at low temperature, m(T) ss 1 — kT/Tq, 
where k = Wk^/Jo for classical spin models, here W 
is the Watson integral. Thus, the ratio between the 
slopes of m(T) for each sublattice is directly related to 
the ratio between the exchange values, Jq, as follows, 
k F6 /k Ni = J 0 Nl /J Fc . It is worth noting that the equilib¬ 
rium magnetization as a function of temperature can be 
fitted to the power law m(T) = (1 — T /Tc) fi which in 
turn gives the low temperature limit m(T) = 1 — kT/Tq- 
And more importantly, it gives a link of the dynamics to 
the equilibrium thermodynamic properties through the 
ratio 


r Ni _ A Fe n m K m 


(17) 


Next, we fit the numerically evaluated m(T) curves to 
the power law mFe(Ni)(T) = (1 — T/Tc) K 1 ’ for T < 
0.5 Tq- This allows us to directly estimate the ratio be¬ 
tween the relaxation times for the three alloys, see Table 
pi} We can see that the relation in Eq. Jl7| agrees well 
for the three alloys even for PyeoCmtg. 

For a more general case, for instance at elevated tem¬ 
peratures, where the one-exponential solution is not a 
good approximation, we have to solve numerically for the 
coefficients of each exponential decay A e and B e . Apart 
from the exchange interactions and temperature depen¬ 
dence, A e and B e also depend on the initial conditions 
<5m e (0) = m e ( 0) — ml. 


E. Effect on distinct local damping parameters on 
the magnetization dynamics 


The intrinsic (atomistic) damping parameters A e are 
not be necessarily the same for both sublattices. To in¬ 
vestigate the effect of different damping parameters we 
consider that the magnetic system is initially at equi¬ 
librium at room temperature T = 300 K. Then a heat 
pulse Tp U i se is applied for 1 ps. We define r Fe(Nl) the 
time at which the normalized magnetization, m norm (t) = 
(m(t) — TO m i n )/(m({ = g) — TO m i n ) is 1/e. The results for a 
broad parameter space of A Fe /A Nl and heat pulse tem¬ 
perature Tp U i se (scaled to Tq) are shown in Fig. 
reffig:PhaseDiagramRelaxTimesPy. The line where 
t Ni /t Fc = 1 lies at low pulse temperature (linear limit 

in the LLB) A Fe /A Nl = 1.563. The critical ratio (yrn) 

is close to the one which could be predicted from Eq. 
fllT) assuming r Nl /r Fe = 1: 

/ A Fe \ n Fe K Fe 

(18) 

Estimations of this critical ratio at low temperatures can 
be found in Table mi The ratio is around 2 for all the 
alloys. 
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_Ni / Fe 
T /T 



FIG. 8. (Color online) Ratio between the relaxation times r 
of the Fe and Ni sublattices in Py after the application of a 
heat pulse of temperature T pu i se for a range of values of the 
ratio of intrinsic damping parameters, A Fc /A Fc . Black lines 
represent A Fo /A Nl values where the ratio between relaxation 
times r Nl and r Fc is constant with the value given by the label. 

The results presented in Fig. [8] show a variety of possi¬ 
ble situations that can be encountered in experiments on 
alloys with two magnetic sublattices. They show that in 
the case of equal coupling to the heat-bath, the Ni sub¬ 
lattice demagnetizes faster than the Fe sublattics in all 
temperature ranges. The situation may be changed if Fe 
is as least twice stronger coupled to the heat-bath than 
Ni. This conclusion is not inconsistent with the dispro- 
portional couplings that were assumed in Ref. [XT] Thus, 
Fe can demagnetize faster than Ni (as reported in Ref. 
IM|l only if Fe is stronger coupled to the heat-bath. 

IV. DISCUSSION AND CONCLUSION 

Element-specific magnetization dynamics in multi- 
subla ttice m agnets has attracted a lot of attention 
lately 124 1 45146 ! The case of GdFeCo ferrimagnetic al¬ 
loys is paradigmatic since this was the first material 
where the so-called ultrafast all-optical switching (AOS) 
of the magnetization has been observecP. The element- 
dependent magnetization dynamics in GdFeC o alloys has 
meanwhile been thoroughly studied 9 10 12 2 ®SD. F rom a 
fundamental view point, however, it is also important to 
understand the element-specific magnetization dynamics 
in multi-element ferromagnetic alloys. This is challenging 
from a modeling perspective and, moreover, contradict¬ 
ing results have been observed in NiFe alloy^SMl. 

To treat such alloys we have developed here a hierar¬ 
chical multiscale approach for disordered multisublattice 
ferromagnets. The electronic structure ab initio calcula¬ 
tions of the exchange integrals between atomic spins in 
FeNi alloys serves as as an accurate foundation to define 


a classical Heisenberg spin Hamiltonian which in turn has 
been used to calculate the element-specific magnetization 
dynamics of atomic spins through computer simulations 
based on the stochastic LLG equation. Our simulations 
predict consistently a faster demagnetization of the Ni 
as compared to the Fe. These findings are however in 
contrast to the dynamics measured by Mathias et al 

From a modeling perspective, we have linked informa¬ 
tion obtained from computer simulations of the atom¬ 
istic Heisenberg Hamiltonian to large scale continuum 
theory on the basis of the recently derived finite temper¬ 
ature LLB model for two sublattice magnetipS. The LLB 
model is rather general, it can be applied not only to fer¬ 
romagnetic alloys, as we have done in the present work, 
but also to ferrimagnetic alloys.^ Thanks to analytical 
expressions coming from the LLB model we have been 
able to interpret the distinct element-specific dynamics 
in FeNi alloys in terms of the strength of the exchange 
interaction acting on each sublattice. Assuming equal 
damping parameters for Fe and Ni, the difference is not 
only coming from the different atomic moments. Analyt¬ 
ical expressions derived for the ratio between demagneti¬ 
zation times in Fe and in Ni compare very well to numer¬ 
ical results from computer simulations of the atomistic 
spin model. To investigate the effect of different intrin¬ 
sic damping parameters we have restrained ourselves to 
use the LLB approach which is computationally less ex¬ 
pensive than the atomistic spin dynamic simulations on 
a large system of atomic spins. Our investigation thus 
prepares a route to an easier characterization, prediction 
and hence, control of the thermal magnetic properties 
of disordered multi-sublattice magnets, something which 
will be valuable for technological purposes. 

As for the applicability of our multiscale approach to 
ferrimagnetic materials, one would obviously need accu¬ 
rately calculated exchange integrals as a starting point. 
Computing these for rare-earth transition metals alloys 
might not straightforward, as the rare-earth ions con¬ 
tain mostly localized /-electrons with a sizable orbital 
contribution to the atomic moment. However it is ex¬ 
pected that for ferrimagnetic alloys, or multilayers with 
antiparallel alignment, composed of transition metals 
this task will be easier. Initial theoretical comparisons 
of the element-specific demagnetization in GdFeCo were 
done recently by Atxitia el al. 10 who obtained a good 
agreement with experimental observations. However, in 
this work the exchange integrals as well as the magnetic 
atomic moments were taken from phenomenological con¬ 
siderations contrary to the present work where all the pa¬ 
rameters are obtained from first-principles calculations. 
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